

### Post-treatment bias simulation
set.seed(14627)
x <- rnorm(500, 50, 15)
z <- rnorm(500, 50, 15)
w <- rnorm(500, 0.5*x + 0.5*z, 5)
y <- rnorm(500, 75+ -0.5*z, 5)

sub <- w > 60 & w < 70


#cairo_pdf(filename = "../figs/post-treatment.pdf", family = "Minion Pro", height = 3.5, width = 6.5, pointsize = 9)
par(mfrow = c(1,2))
plot(x, z, type = "n", pch = 19,
     xlab = "Treatment", ylab = "Confounder", yaxt = "n", bty = "n",
     xlim = c(0,100), ylim = c(0,100), las = 1)
axis(side = 2, las = 2)
points(x[!sub], z[!sub], col = "grey60", cex = 0.75)
abline(lm( z[!sub] ~ x[!sub]), col = "grey60")
points(x[sub], z[sub], pch = 19, col = "black", cex = 0.75)
abline(lm( z[sub] ~ x[sub]))
text(x = 50, y = 100, expression("mediator" %in% group("(", list(60, 70), ")")))
text(x = 50, y = 2.5, expression("mediator" %notin% group("(", list(60, 70), ")")), col = "grey60")
plot(x, y, type = "n", pch = 19,
     xlab = "Treatment", ylab = "Outcome", yaxt = "n", bty = "n",
     xlim = c(0,100), ylim = c(0,100), las = 1)
axis(side = 2, las = 2)
points(x[!sub], y[!sub], col = "grey60", cex = 0.75)
abline(lm( y[!sub] ~ x[!sub]), col = "grey60")
points(x[sub], y[sub], pch = 19, col = "black", cex = 0.75)
abline(lm( y[sub] ~ x[sub]))
text(x = 60, y = 15, expression("mediator" %in% group("(", list(60, 70), ")")))
text(x = 45, y = 80, expression("mediator" %notin% group("(", list(60, 70), ")")), col = "grey60")
#dev.off()

nsims <- 10000
truth <- rep(NA, times = nsims)
biased <- rep(NA, times = nsims)
seqg <- rep(NA, times = nsims)
for (i in 1:nsims) {
  x <- rnorm(500, 50, 15)
  z <- rnorm(500, 50, 15)
  w <- rnorm(500, 0.5*x + 0.5*z, 5)
  y <- rnorm(500, 75+ -0.5*z, 5)

  sub <- w > 60 & w < 70
  truth[i] <- coef(lm(y ~ x))["x"]
  biased[i] <- coef(lm(y ~ x + w))["x"]

  first <- lm(y ~ x + w + z)
  seqg[i] <- coef(lm(I(y - coef(first)["w"]*w) ~ x))["x"]
}

#cairo_pdf(filename = "../figs/ptbias-sims.pdf", family = "Minion Pro", height = 3, width = 4.5, pointsize = 9)
plot(density(truth), col = "grey60", xlim = c(-0.2, 0.5), ylim = c(0,20), lwd = 2, bty = "n",
     yaxt = "n", xlab = "Coefficient on treatment", main = "")
axis(side = 2, las = 2)
abline(v = 0, lty = 2, col = "grey60")
lines(density(biased), col = "black", lwd = 2)
lines(density(seqg), col = "dodgerblue", lwd = 2)
text(x = 0, y = 20, "True effect", pos = 4)
text(x = 0.375, y = 10, "Conditioning\n on M", pos = 4)
text(x = 0.025, y = 10, "Unconditional", pos = 4, col = "grey60")
text(x = 0.05, y = 5, "Sequential\ng-estimation", pos = 4, col = "dodgerblue")
#dev.off()
